Motivação do classificador de câncer de mama

O objetivo desse notebook é aplicar técnicas de classificação aos dados do dataset Breast cancer Wisconsin, disponível aqui. O interesse nesse dataset vem do outubro rosa, mês de conscientização para informar, alertar e previnir o câncer de mama. De acordo com o Instituto Nacional de Câncer (INCA), em 2019 foram estimados 59700 novos casos da doença no Brasil.

Quanto antes, melhor!

Esse estudo é apenas para fins didáticos e de aprendizado, que não necessariamente reflete a realidade ou uma opinião médica.

Pacotes

library(tidyverse)
library(janitor)
library(GGally)
library(caret)
library(printr)
library(corrplot)
library(factoextra)
library(pROC)

set.seed(889)

source("analysis.r")

## [1] "AUC on Test 0.999689440993789"

Explorando os dados: Correlações e separabilidade das classes

Mean

É possível ver rapidamente que algumas variáveis estão correlacionadas fortemente entre si.

Algumas variáveis apresentam correlação positiva quase perfeita (muito próximo de 1), como

  • radius_mean e perimeter_mean (0.998)

  • radius_mean e area_mean (0.987)

  • perimeter_area e area_mean (0.987)

Outras variáveis também apresentam uma correlação forte (acima de 0.5), como

  • concavity_mean e compactness_mean (0.883)

  • concave.points_mean e concavity_mean (0.921)

  • concave.points_mean e perimeter_mean (0.851)

Quanto à separabilidade dos dados, observando os histogramas, parece que as variáveis radius_mean, perimeter_mean, area_mean e concavity_mean ajudam a melhor identificar as populações benigno e maligno, mas ainda assim vemos sobreposições e nenhuma variável parece separar perfeitamente essas classes.

Standard Error

As correlações do erro padrão dessas variáveis são, de maneira geral, um pouco mais baixas, mas ainda sim temos correlações fortes nas mesmas variáveis da média.

Quanto à separabilidade das classes, as variáveis parecem não ajudar a identificar as duas populações.

Worst

Os scatterplots dessas variáveis, que são as médias dos maiores valores, também são bem parecidas com as observadas para as médias.

Por meio das variáveis concave.points_worst, radius_worst e perimeter_worst parecem ser as que mais ajudam separar as classes do diagnóstico.


Quando uma variável tem correlação perfeita, sabendo o valor da variável X, conseguimos prever o de Y. Com isso, podemos escrever X em função de Y (ou vice-versa). E observamos nas matrizes acima que algumas variáveis têm correlação bem alta.

Trabalhar com dados altamente correlacionados:

Para ajudar nesse problema, vamos recorrer à técnicas de redução de dimensionalidade.

PCA: Reduzindo a dimensão

Análise de componentes principais é uma técnica de fatorização de matriz em que é possível explicar (parte da) a variabilidade dos dados por meio de combinações lineares não correlacionadas dos dados originais. De maneira geral, a ideia de fatorar uma matriz é escrever uma matriz “complicada” no produto de duas matrizes “mais simples”.

Vamos definir as componentes principais por meio de operações matriciais, para conhecer todo o processo.

  1. Estimando a matriz de correlação:

O primeiro passo é estimar a matriz de correlação dos nossos dados. Na verdade, nessa etapa usa-se a matriz de covariâncias, mas um problema que pode surgir é que algumas variância são maiores que outras, e em casos muito discrepantes isso distorce as componentes principais. Para evitar isso, o procedimento é: padronizar a matriz de covariâncias e estimar as componentes a partir dessa matriz padronizada. Mas isso é o mesmo que usar a matriz de correlação e por isso vou seguir por ela.

upper <- cor_matrix
upper[upper.tri(cor_matrix)] <- NA 

knitr::kable(upper[1:5, 1:5],  caption = "Matriz de correlação: 5 primeiras linhas/colunas")
Matriz de correlação: 5 primeiras linhas/colunas
radius_mean texture_mean perimeter_mean area_mean smoothness_mean
radius_mean 1.0000000 NA NA NA NA
texture_mean 0.3237819 1.0000000 NA NA NA
perimeter_mean 0.9978553 0.3295331 1.0000000 NA NA
area_mean 0.9873572 0.3210857 0.9865068 1.0000000 NA
smoothness_mean 0.1705812 -0.0233885 0.2072782 0.1770284 1
  1. Obtendo os autovalores:
eig$values
##  [1] 1.328161e+01 5.691355e+00 2.817949e+00 1.980640e+00 1.648731e+00
##  [6] 1.207357e+00 6.752201e-01 4.766171e-01 4.168948e-01 3.506935e-01
## [11] 2.939157e-01 2.611614e-01 2.413575e-01 1.570097e-01 9.413497e-02
## [16] 7.986280e-02 5.939904e-02 5.261878e-02 4.947759e-02 3.115940e-02
## [21] 2.997289e-02 2.743940e-02 2.434084e-02 1.805501e-02 1.548127e-02
## [26] 8.177640e-03 6.900464e-03 1.589338e-03 7.488031e-04 1.330448e-04

plot(eig$values)
abline(h = 1)

Para estimar o número de componentes principais ideal, um critério prático é observar os autovalores maiores ou iguais a 1, pois dessa forma as combinações lineares explicam a variância de uma variável original (padronizada) do dataset. Observando o scree-plot acima, é k = 6 é o número de componentes principais adequados para explicar a estrutura da variância dos dados.

  1. Obtendo os autovetores e escrevendo a 1ª componente principal:

Observando os autovetores associados aos seis autovalores identificados (apenas as 6 primeiras linhas de 30)

V1 V2 V3 V4 V5 V6
-0.2189024 -0.2338571 -0.0085312 0.0414090 0.0377864 0.0187408
-0.1037246 -0.0597061 0.0645499 -0.6030500 -0.0494689 -0.0321788
-0.2275373 -0.2151814 -0.0093142 0.0419831 0.0373747 0.0173084
-0.2209950 -0.2310767 0.0286995 0.0534338 0.0103313 -0.0018877
-0.1425897 0.1861130 -0.1042919 0.1593828 -0.3650885 -0.2863745
-0.2392854 0.1518916 -0.0740916 0.0317946 0.0117040 -0.0141309

Olhando um pouco para a teoria, a j-ésima componente principal é dada por:

\[ Ŷ_j = ê_{j1}X_{1} + ê_{j2}X_{2} + ... + ê_{jp}X_{p} \]

E podemos escrever, para os dados do dataset em questão, a primeira componente principal da seguinte forma:

\[ Ŷ_1 = -0.2189\text{radius_mean} -0.1037\text{texture_mean} - 0.227\text{perimeter_mean} + ... - 0.132\text{fractal_dimension_worst} \]

E assim em diante…

Mas como visualizações são mais interessantes, vamos construir gráficos!

Visualizando as componentes principais

Bom, como a função fviz_pca_biplot do pacote factoextra faz um gráfico apropriado a partir de um objeto da classe PCA, vou utilizar a função princomp (do stats, que faz parte do R base) para criar esse objeto. Essa função foi escolhida porque ela usa o teorema da decomposição espectral, que é a abordagem teórica utilizada na descrição acima. Caso a abordagem fosse decomposição em valores singulares (SVD), a função apropriada seria a prcomp (também do stats).

pca <- princomp(data_processed[,-1], cor = TRUE, scores = TRUE)

fviz_pca_biplot(pca, col.ind = data_processed$diagnosis, col="black",
                palette = "Dark2", geom = "point", repel=TRUE,
                legend.title="Diagnóstico", addEllipses = TRUE) +
  theme(axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) 

A primeira componente explica 44.3% da variância dos dados. É possível observar que existem variáveis correlacionadas negativamente, como smoothness_se e radius_mean, justamente os vetores mais distantes entre si. Também vemos que há ali um agrupamento das classes Benigno/Maligno.

Como o PCA foi utilizado para reduzir a dimensionalidade dos dados, os scores das 6 primeiras dimensões serão utilizados para ajustar o modelo de classificação.

Classificando

Organizando os dados que serão utilizados no modelo:

data_scores <- 
  bind_cols(data_processed$diagnosis, as.data.frame(pca$scores[,1:6])) %>% 
  clean_names() %>% 
  rename(diagnosis = x1)
## New names:
## * NA -> ...1

Dividindo os dados em conjuntos de treino e teste (70/30):

train <- data_scores %>% sample_frac(.7)
test <- data_scores %>% anti_join(train)

Modelo 01: Regressão logística

Ajustando o modelo

logistic_model1 <- glm(diagnosis ~ ., family = binomial, data = train)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(logistic_model1)
## 
## Call:
## glm(formula = diagnosis ~ ., family = binomial, data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6845  -0.0496  -0.0045   0.0012   3.4221  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.2843     0.3626  -0.784  0.43300    
## comp_1        2.9268     0.5107   5.731 9.96e-09 ***
## comp_2        1.6343     0.3493   4.679 2.89e-06 ***
## comp_3        0.4722     0.2297   2.056  0.03982 *  
## comp_4       -0.7846     0.2670  -2.939  0.00329 ** 
## comp_5       -1.3397     0.4449  -3.011  0.00260 ** 
## comp_6       -0.3746     0.2821  -1.328  0.18411    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 533.015  on 397  degrees of freedom
## Residual deviance:  57.436  on 391  degrees of freedom
## AIC: 71.436
## 
## Number of Fisher Scoring iterations: 9

A variável comp_6 parece não contribuir significativamente para classificarmos o diagnóstico (utilizando \(\alpha = 0.05\)). Um novo modelo pode ser ajustado removendo essa variável, já que variáveis irrelevantes podem deteriorar a taxa de erro do modelo.

logistic_model2 <- glm(diagnosis ~ comp_1 + comp_2 + comp_3 + comp_4 + comp_5, 
                       family = binomial, data = train)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(logistic_model2)
## 
## Call:
## glm(formula = diagnosis ~ comp_1 + comp_2 + comp_3 + comp_4 + 
##     comp_5, family = binomial, data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7274  -0.0542  -0.0060   0.0016   3.4014  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.2408     0.3643  -0.661  0.50861    
## comp_1        2.8145     0.4734   5.945 2.77e-09 ***
## comp_2        1.5617     0.3314   4.713 2.44e-06 ***
## comp_3        0.3441     0.1971   1.746  0.08086 .  
## comp_4       -0.7896     0.2631  -3.002  0.00269 ** 
## comp_5       -1.2132     0.3975  -3.052  0.00227 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 533.015  on 397  degrees of freedom
## Residual deviance:  59.161  on 392  degrees of freedom
## AIC: 71.161
## 
## Number of Fisher Scoring iterations: 9

E como fica no conjunto de teste?

Queremos prever novas observações, então vamos lá:

pred <-
  predict(logistic_model2, test, type = "response") %>% 
  as.data.frame() %>% 
  rename(., "prob" = ".") %>% 
  bind_cols(test) %>% 
  mutate(predicted = as.factor(ifelse(prob > 0.5, "Malignant", "Benign")))

As probabilidades são dadas na forma \(P(Y = 1 | X)\), e olhando para a dummy definida como 1

contrasts(data_scores$diagnosis)
Malignant
Benign 0
Malignant 1

olhamos para os valores como sendo a probabilidade do tumor ser maligno.

Avaliando a acurácia do modelo

cmat <- conf_mat(table(pred$predicted, pred$diagnosis))
autoplot(cmat, type = "heatmap")


pred %>% metrics(diagnosis, predicted)
.metric .estimator .estimate
accuracy binary 0.9766082
kap binary 0.9473765

A acurácia do modelo é de 99.4% e o coeficiente de Kappa 98.7%, que nos indica uma concordância entre o predito e os valores verdadeiros. É importante lembrar que a acurácia é uma métrica que deve ser usada com cuidado, pois ela pode não representar bem a situação.

Na área da saúde, as medidas mais utilizadas são a especificidade e a sensibilidade, comuns em testes diagnósticos. Em termos práticos, a sensibilidade avalia a capacidade do classificador detectar o tumor maligno quando ele está presente. Já a especificidade nos dá a probabilidade de um teste dar negativo (no caso, benigno) quando, de fato, não há doença. A sensibilidade do classificador é de 98,2% e a especifidade 100%.

Quando estamos em um cenário de classificação, uma forma de checarmos a taxa de erro associada ao conjunto de teste é estimada da seguinte forma:

\(Media(I(y_{0} \neq \hat{y_0}))\)

Ou seja, é a média de uma variável indicadora que assume 1 quando a classe predita é igual a classe verdadeira e 0, caso contrário. Queremos o menor valor possível para essa média. Assim, considerando o problema, temos:

mean(pred$diagnosis != pred$predicted)
## [1] 0.02339181

Ou seja, o modelo de regressão logística previu incorretamente, em média, apenas 0.6% das vezes.

ROC e AUC:

predr <- prediction(pred$prob, pred$diagnosis)
perf <- performance(predr, "tpr", "fpr")
plot(perf)



auc <- performance(predr, measure = 'auc') 
auc <- unlist(slot(auc, 'y.values'))
print(paste('AUC on Test', auc))
## [1] "AUC on Test 0.996428571428571"

Referências

https://www.inca.gov.br/publicacoes/livros/situacao-do-cancer-de-mama-no-brasil-sintese-de-dados-dos-sistemas-de-informacao